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^ '. Abstract 

Q, 

1^^ , The transient thermocapinary migration of drops with nontrivial deformation is studied. The 

finite difference method is employed to solve the incompressible Navier-Stokes equations coupled 

J>^. with the energy equation; the front-tracking method is adopted to track the moving deformable 

^ ■ 
^ i drop interface. In the hot region, deformations of drops increase with the decrease of interfacial 

tensions. In order to indicate the temperature impact on the interfacial tension, a local capillary 

number (Cai) is introduced. It is found that, when the drop density is smaller/larger than that of 

the bulk fluid, the drop velocity decreases/increases with the increase of the drop deformation. 
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I. INTRODUCTION 



Thermocapillary migration, which is induced by the variance of the interfacial tension 
when a temperature gradient is imposed on the bulk fluid, is an important phenomenon in 
drop/bubble transportation. With the development of aeronautics and astronautics, it is 
necessary to study the thermocapillary motion because of its practical roles in the material 
processing, and the management of heat and fluids in space [l|, y]- In most studies on the 
thermocapillary motion, it is assumed that the interfacial tension linearly depends on the 
temperature: 

a(r) = oro + aT(r-To), (1) 

where the constant o"t is the interfacial tension coefficient of temperature, and o"o the inter- 
facial tension at the reference temperature Tq. 

Considering the balance between the capillary force and viscosity on the bubble or drop, 
the reference velocity can be defined as 



U = |orr||V Too|-Ro//^i, 

in which |V Tod is the temperature gradient imposed on bulk fiuid (the subscript 1 repre- 
sented), /i the viscosity, and Rq the radius of the spherical drop/bubble. 





ao{W-^N/m) 


aT{W^N/{m • K)) 


Gas bubble/Silicone oil[3] 


17.3 


-0.061 


Fc-75 drop/Silicone oil [3] 


3.47 


-0.036 


Monotectic alloy [4] 


5.7 


-0.23 



TABLE I: Typical interfacial tensions and their temperature coefficients. 



The capillary number is defined as 



Ca = fiiU/(7o = IcttIIVTooI-Ro/o^o- 

In most cases, the Ca numbers in thermocapillary migration are small (0.01 or less, e.g., 
the materials in the first two lines of Table 1), and the deformations of drops are also 

n 

very small. Therefore, previous studies usually assumed drops to be nondeformable [5[. 
However, in practice, there are many situations where the interface tension is very small, 



or where the deformation introduced by capillary effect is very large. For the case of the 
preparation of monotectic alloys in micro-gravity conditions [6|, the Ca number of the Al- 
based alloy {Al34^Bif^^^^)Q^Sn^ is 0.12 at the temperature of 1200k if Rq = 300/im and 
|VToo| = lOOK/cm [7] (see the third line in Table 1). Such a big Ca number can cause 
obvious deformation on the drop. Similar situations happen in the chemical flooding methods 



applied in enhanced oil recovery |8l-ll2|. Finally, experiments in space also observed some 
slight deformation of fairly large bubbles [3]. Therefore, it is necessary to investigate the 
deformed drops with relatively large Ca's or small interfacial tensions in the thermocapillary 
studies. 

When inertia and thermal convections are neglected, the velocity of the spherical drop in 
thermocapillary motion is presented in the pioneering work of Young et al. jlSj : 

Vygb = 2AU, 

where 

"^ " (2 + 3a)(2 + A)' ^^^ 

and a = yU2//^i and A = ^2/^1 are the ratios of the kinematic viscosity and thermal conduc- 
tivity between droplet (the subscript 2 represented) and background liquid, respectively. 
There are a lot of investigations studying the thermocapillary motion, which are reviewed 
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15| . Most of the researches sofar focus on the convection effects of inertia and energy, 
and only a few of them are about the drop/bubble deformation. When the thermal convec- 
tion is very small, the exact solution of the momentum equation was provided in [16|, where 
small inertial deformations of drops were also calculated. Under the same assumption, using 
the Lorentz reciprocal theorem, Haj-Hariri et al. calculated the small inertial deformations 
as well as the consequent correction on the temperature field and the migration velocity 



17(1 . Later, the work is extended to moderate parameters in a steady-state numerical inves- 



tigation 



1^ 



It is concluded that even a small deformation will retard the motion of the 
drop, and the steady-state migration can still be reached even for large capillary numbers 
up to 0.5. However, the actual continuous decrease of the interfacial tension when the drop 
migrates to the hot region was not considered in their study. While taking this decrease into 
consideration, bubbles with large capillary numbers can not reach steady migrating states 
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FIG. 1: The sketch of thermocapinary motion m the axisymmetric model. 

According to our knowledge, transient migration behaviors due to drop deformation have 
never been reported before, which is the main topic in this work. This paper is arranged 
as follows: the governing equations and numerical methods are introduced in Section 2, the 
validation tests are presented in Section 3, and the results and discussions are in Section 4. 



II. GOVERNING EQUATIONS AND NUMERICAL METHODS 

In our axisymmetric simulation, the initially-spherical droplet with a radius of Rq is 
surrounded by the bulk liquid, which is imposed on a constant temperature gradient along 
the z axis (Fig. 1). The calculated domain is Q{{r,z)) (r G [0,ri], z G [0,2i]), and the 
center of the droplet is at r = 0. The governing equations for the thermocapillary migration 
of droplets are: 

V ■ u = 0, (3) 



dt 



+ V ■ (puu 
,dT 



-Vp + V ■ (/i(Vu + V' u)) + F, 



pCp(— + u-VT)=V-(A;VT), 



(4) 
(5) 



in which u = {vr,Vz) is the velocity, p the density, p the pressure , and Cp the specific heat. 



F^ is the body_ 



presented in [20 



' orce term due to the interfacial tension, and its axisymmetrical formula is 

da 



2l|: 



/ (5(x — Xf)(cr/tn + — -T)(is. 

JB OS 



(6) 



Here, x = (r, z) is the space vector, x/ = (r/, Zf) the position of the cell / on the interface 
B, and s the natural coordinate along the interface, n = {nr,nz) and r denote the normal 
and tangential unit vectors of the interface, respectively, k = ki + K2 is the sum of two 
principal curvatures of the interface, Ki is the in-plane curvature and ^2 is given as [22| : 



K2 = • 

We defined nondimensional variables as follows: 



(7) 



,i?0^ 



u = u/U, X = x/i?o, t = t/{—) 



f = T/(|VToo|i?o), P = p/pu ^ = p2/pi, p = p/pi, (8) 

k = k/ki, Cp = Cp/Cpi, 7 = Cp2/Cpi, 
F, = F.i?o/(pif/'), Ma = p,CpiURo/k,. 
In total, the problem is governed by seven nondimensional parameters: 



Re, Ma, Ca, a. A, 7 and ^. The nondimensional equations in the axisymmetrical model are: 

0, (9) 
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(12) 



Using Eq. [H the nondimensional interfacial tension term F„ can be written as: 



(13) 



Eqs. [9IIT2] are valid for both phases of the drop and the bulk fluid. The discontinu- 
ities of the physical parameters across the interface are handled with an indicator function 
for more details). The governing equations are discretized by the second-order 



see 



center- difference method on fixed, staggered, Cartesian grids. The projection method 24 1 
is adopted to solve the Navier-Stokes equations. 

On the symmetric axis, the symmetrical condition is employed: 

dvz df 

Vr\f=0 = 0, — ^|,-.=o = 0, —^1^=0 = 0, (14) 

or or 

and all other boundaries for velocities adopt rigid boundary conditions, and those for tem- 
perature use constant temperature boundary conditions. Following the tradition in this 
field, the initial conditions are: 

t;r|f=o = ^^|f=o = 0, f\t=o = z. (15) 

In the following, symbols without bars will be used to denote nondimensional values. 

In addition, when the momentum equation is solved at each new time step, the front 
element / is supposed to have the same velocity (uj) as that of the surrounding fluid. In 
order to update the front with a conserved drop volume, a velocity correction Am„ is imposed 
on Uf in the normal direction (n/) of the element: 

With this correction, the drop volume is almost conserved throughout our simulations with 
a maximum volume loss of less than 0.1% of the initial volume. By comparison, without the 
velocity correction, there will be about 2% loss at the simulating time t = 20. 

The computing domain in this paper is {r,z) G [0,6] x [0,24], and the resolution is 
128 X 512. The time step is 2 x 10~^ for most cases, and 1 x 10~^ for the ^ = 0.01 case in 
subsection 4.2. 

III. THE VALIDATION OF THE COMPUTING PROGRAM 

In tradition, the scaled deviation of the drop profile from the sphere is defined as 

f{9) = R{9)/Ro - 1, 

where, the polar angle 6 is measured from the front stagnation point and R{0) denotes the 
distance between the interface and the mass center of the drop (Fig. 1). 
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FIG. 2: Deviation of the drop profile from the sphere for Ma = 1, Ca = 0.1, a = A = ^ = 0.5, 

[171]; sohd hne: 



7 = 0.25, and Re as indicated. Dash line: the theoretical prediction in 
work. 
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FIG. 3: Deviations of drop profiles from spheres in [13], [iSj], and this work, where Re 
Ma = l,Ca = 0.1, a = X = ^ = 0.5, and 7 = 0.25. 
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Before this work, there was a theoretical result of small inertial deformations of drops 
with small Ca and Re numbers [l7|: 

f(9) = lA^ReCa{^-l){3cos^9-l). (17) 

however, used a different formula for the interfacial 



The earlier numerical simulation 
tension: 



Compared with Eq. [13], the effect of capillary force in Eq. [18] was assumed constant on the 
normal direction of the interface, or, there is no decrease in the interfacial tension when the 



droplet moves to the hotter region 18|. There are two distinguished differences between 
these two assumptions: 

1. As it will be shown later in this paper, Eq. [13] leads to a bigger drop deformation than 
Eq. [13 

2. With Eq. [T31 the simulation must be stopped before the drop moves a distance of 
(1/Ca — 1) in the z direction. This is because the interfacial tension of the drop will 
be smaller than zero beyond that location, and not only the simulation will collapse, 
but also it is not physically possible. 

With Eq. [181 on the other hand, the simulation can be extended to the infinity. 

In this section, Eq. [18] is adopted in our codes to have some comparisons with the 
previous researches. It seems that our simulations have very good agreements with the 
analytical work (Fig. 2), and that a better similarity is achieved when the Reynolds number 



is smaller. Compared with the previous numerical work 18|, our numerical simulations are 
closer to the theoretical analysis (Fig. 3). 

In the rest of this paper, the interfacial tension is calculated by Eq. [131 cind simulations 
are stopped before the capillary force on the drop becomes negative. 

IV. RESULTS AND DISCUSSIONS 

A. Transient thermocapillary migrations of droplets with large capillary numbers 

In this subsection, we simulate the thermocapillary process with the same set of param- 
eters in |18j: Re = Ma = 50, Ca = 0.1, a = X = C, = 0.5, and 7 = 0.25. The time evolution 
of migration velocity is shown in Fig. 4 (the solid line). When the capillary number is suffi- 
ciently small {Ca = 0.05, the dashed line in Fig. 4), the drop will reach a steady migrating 
state after the initial accelerating process. However, the velocity with Ca = 0.1 has an 
accelerating decrease when t > 60, and can not reach a steady state. This result is different 
in which the steady velocity state was assumed in advance. 



from that of 



The deviation of the drop profile from the sphere for Ca = 0.1 at different times is plotted 
in Fig. 5. It seems that when the drop migrates to the hotter region, the deformation keeps 
increasing and the drop loses its fore and aft symmetry. At the late stage, the bottom of 
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FIG. 4: Time evolution of drop velocity for FIG. 5: Deviation of the drop profile from the 
i?e = Ma = 50, a = A = ^ = 0.5, and 7 = 0.25. sphere at indicated times for Re = Ma = 50, 
Ca = 0.1 (solid line); Ca = 0.05 (dashed line). Ca = 0.1, a = A = ^ = 0.5, and 7 = 0.25. 

the drop becomes flattened, and the drop looks like a cap (Fig. 6(b)). The deformation 



18|, 



at t = 60 in the current study is about twice bigger than that in the simulation of 
and even bigger difference at the late stage. There are two reasons for the rapid decline in 
velocity: 

1. Because the drop becomes flatter, the resistance on the drop is increased. 

2. Because the distance between the front and rear stagnation points at the late stage is 
shorter than that in the beginning, their temperature difference also decreases (Fig. 
7). As a result, the thermocapillary driving force on the drop becomes smaller. 

We also studied the influence of different Re numbers. Ma numbers, and Ca numbers. 
Some typical results are shown in Figs. 8. With other parameters fixed, the migrating 
velocity tails off earlier for the larger Reynolds number. 

With Re = Ma = 50, Ca = 0.1, and ^ = 0.5, we studied the influence of a, A and 7. 
These three parameters play trivial roles in the transient motions of deformable drops. With 
small a, A or 7, drops migrate faster (see p, |25|) and need shorter time to reach the hot 
regions with near- zero interfacial tensions. As a result, the migration velocities tail off at 
earlier times. 

To sum up, the influence of all other parameters except ^ depends on when the drop 
migrates to the warm region with very low interfacial tension, in other words, the speed of 
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FIG. 6: Isotherms around the drop for Re = Ma = 50, Ca = 0.1, a = X = ^ = 0.5, and 7 = 0.25. 
a) t = 40; b) t = 70. 
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FIG. 7: The temperature distribution on the drop interface for Re = Ma = 50, Ca = 0.1, 
a = A = ^ = 0.5, and 7 = 0.25 when t = 40 and t = 70. 
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FIG. 8: Time evolutions of the drop velocity for different Ca numbers. Ma = 100, a = X = ^ = 0.5, 
and 7 = 0.25. (a) Re = 1; (b) Re = 50. 

the drop. The influence of various parameters on drop speeds has been discussed in detail 

nn 

in our earher work la, um ■ 



B. Influence of density ratio ^ 

In the last subsection, densities of drops are always smaller than those of bulk fluids, and 
we will adopt denser drops {C, = 2) in the following. Several Ca numbers are studied, all 
other parameters are flxed to make the discussion simpler [Re = Ma = 50, a = X = 0.5, 
and 7 = 0.25). The time evolutions of velocities are plotted in Fig. 9, and the most obvious 
difference from Figs. 4&8 is the continuous increasing trend of velocities. To explain this 
phenomenon, we focus on the case of Ca = 0.1: 

1. At the late stage of the migration, the drop proflle is elongated in the z direction (Fig. 
10(b)), and the temperature difference between the front and rear stagnation points 
of the drop becomes larger (Fig. 11). As a result, the thermocapillary driving force 
on the drop becomes larger. 

2. Because the drop becomes slender, the resistance on the drop decreases. 

To further discuss drop deformations, it is essential to introduce the local capillary number 
Cat to indicate the local magnitude of interfacial tension, which is deflned on the moving 
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FIG. 9: Time evolutions of drop velocities for different capillary numbers. Re = 50, Ma = 50, 
a = A = 0.5, 7 = 0.25, and ^ = 2. 
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FIG. 10: Isotherms around the drop at t 
a = X = 0.5, 7 = 0.25, and ^ = 2. 



40 and t = 80. Re = 50, Ma = 50, Ca = 0.1, 
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FIG. 11: Temperature distribution on the drop interface at t = 40 and t = 80. Re = 50, Ma = 50, 
Ca = 0.1, a = X = 0.5, 7 = 0.25, and £ = 2. 
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FIG. 12: (a) Time evolutions of drop velocities with different £'s, the circles on the velocity curves 
indicate the locations at Cai = 0.5; (b) deviations of drop profiles from spheres when Cai = 0.5. 
Re = Ma = 50, Ca = 0.1, a = \ = 0.5, and 7 = 0.25. 
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Ca 



drop center {zc): 

1 - Ca(Tco - To) 1 - Ca{zc - z^ ' 
where, T^a represents the initially temperature at z^.. Time evolutions of drop velocities for 



(19) 



different ^'s are plotted in Fig. 12(a). To investigate the ^ influence on deformation, we 
fixed Cai at the value of 0.5, when the drop starts to have an obviously different velocity 
from that of the non-deformable drop. In the current case [Ca = 0.1), Cai = 0.5 means 
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FIG. 13: Time evolutions of drop profiles for different ^'s, where Re = Ma = 50, Ca = 0.1, 
a = X = 0.5, and 7 = 0.25. From the bottom to the top, the profiles in each column are 
corresponding to the time at t = 0, 20, 40, 60, 80, 100, respectively. 

that the drop migrates a distance of 8 (represented by circles on the curves of Fig. 12(a)). It 
is found that the drop profile deforms to oblate/slender when ^ is smaller/larger than unit. 
When the absolute value of ^ — 1 becomes larger, the deformation is also larger (Fig. 12(b)). 
This is consistent in the analytical result (Eq. [T7|) . except that we have larger deformation 
and RehMa numbers here. Fig. 13 presents shape evolutions of drops with different ^'s. 

It should be noticed that only slight deformations on large bubbles (^ ~ 10^'^) have been 
observed in previous space experiments, and no measurable deformations for those heavy 
drops (^ = 1.98) [Sj. The reason for the difference between experiments and simulations 
might be: 

1. Bubbles/Drops in experiments do not migrate to areas with high temperatures, so 
their interface tensions are still large enough to hold perfect spherical shapes; 

2. The assumption for the cr(T)'s linear dependence on temperature is not correct when 
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FIG. 14: Sketch of normal stress balance on the stagnation point of the drop. 

cr{T) — )■ in experiments (this is, however, the general assumption in most simulations 
in this field). 

To explain the relation between different ^'s and drop-deforming orientations, a simplified 
physical explanation is introduced in the following. Only the front stagnation point is 
studied. 

As shown in Fig. 14, the total pressures are p2 + |P2^^ inside the drop and pi + \piV'^ 
outside the drop; here V is the local fluid velocity, and p2 and pi denote the static pressures 
at the two sides of the stagnation point, respectively. The extra pressure caused by the 
interfacial tension is 2(jH (H is the average curvature). When the interfacial tension is large 
enough, the normal stress balance on the interface is dominated by the interfacial tension 
and the droplet keeps spherical. On the other hand, when the interfacial tension becomes 
small at the hot region (or, Cai > 0.5), the dynamic pressure becomes dominant: 

if ,^ > 1, the dynamic pressure inside the drop (|p2^^) is larger than that outside {^piV'^). 
Larger extra pressure {2aH) is needed to balance the normal stress on the interface, 
so the drop is elongated in the axis direction to have a larger H. 

if ,^ < 1, the dynamic pressure inside is lower than that outside. Smaller extra pressure is 
needed to balance the normal stress, so the drop is compacted in the axis direction to 
have a smaller H. 

V. CONCLUSIONS 

In this work, we found that the influence of deformations on drops is much more com- 
plicated than that on bubbles. When the drop density is smaller than the bulk fluid, the 

15 



migration velocity will decrease with the increase of the deformation. When the drop density 
is larger than the bulk fluid, the migration velocity will increase with the increase of the 
deformation. With the assumption adopted in this paper, we believe that in the hot region 
the drop will always start to deform, and there is no way to have any constant migrating 
velocity. Hence, keeping Ca and Cai small throughout simulations is the only possible way 
for the drop to keep spherical and reach the steady migrating state. 

All the discussions above are based on the assumption of the linear temperature depen- 
dence of the interfacial tension. More complicated temperature dependence will be used 
in the future work, and variations of physical properties with the temperature will also be 
considered. 
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